Study of the island morphology at the early stages of Fe/Mo(110) MBE growth 
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We present theoretical study of morphology of Fe islands grown at Mo(llO) surface in sub- 
monolayer MBE mode. We utilize atomistic SOS model with bond counting, and interactions 
of Fe adatom up to third nearest neighbors. We performed KMC simulations for different values of 
adatom interactions and varying temperatures. We have found that, while for the low temperature 
islands are fat fractals, for the temperature 500A' islands have faceted rhombic-like shape. For the 
higher temperature, islands acquire a rounded shape. In order to evaluated qualitatively morpho- 
logical changes, we measured averaged aspect ration of islands. We calculated dependence of the 
average aspect ratio on the temperature, and on the strength of interactions of an adatom with 
neighbors. 

PACS numbers: 68.35.Fx, 68.55.A-, 81.15.Aa 
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I. INTRODUCTION 

Iron films grown on molybdenum are known to have in- 
teresting ferromagnetic properties [ij. Hence, nature and 
formation of the iron islands at the early stages of molec- 
ular beam epitaxy (MBE) growth are attracting great 
amount of researchers' attention. Knowing the processes 
of island formation in sub-monolayer regime, will provide 
film manufacturers the ability to control the growth pro- 
cess, predict the outcome of deposition and produce films 
with needed properties. 

A number of experimental investigations were done on 
this subject. For the review of experimental works on the 
morphology and magnetic behavior of the iron islands we 
address the reader to [ll, [1] and references therein. 

Theoretical aspect of the early stages of MBE growth 
on bcc(llO) surface still remains an open question. Some 
work has been done to develop a model that repro- 
duces some experimental results for homoepitaxy on 
Fe(llO) 3. Recently, results of the phase-field simula- 
tion of stripe arrays growth on bcc(llO) surface in sub- 
monolayer regime has been reported But, to our 
knowledge, there are no deep theoretical investigation of 
the problem of island formation for Fe/Mo(110) system. 

The main goal of this work is to study morphology 
of the islands in the sub-monolayer regime. Here, we 
present some preliminary results on the island forma- 
tion at early stages of Fe/Mo(110) MBE growth. The 
structure of the paper is as follows. In Section III Al we 
give short theoretical description of MBE growth model. 
To implement the model we are using on-lattice Kinetic 
Monte Carlo (KMC) algorithm, which is widely used for 
modelling of MBE growth. In the beginning of Sec- 
tion nil Al we define quantities utilized for characteri- 
zation of island morphology. Then this section contin- 
ues by description of changes of island morphology while 
varying substrate temperature. Quantitative parameter 



(average aspect ratio) as well as concept of island com- 
pactness are used. In Section IIII Bl we study the mor- 
phology of the island as a function of different values of 
interaction energies. We present results for the depen- 
dence of island aspect ratio on interaction energies with 
nearest and second nearest neighbors. Each section also 
contains visual illustrations of islands shape for different 
parameters. Section llVl concludes our current results and 
outlines the future directions of work. 



II. MODEL 
A. Growth Model 

In our study, we used growth model based on solid-on- 
solid (SOS) model of molecular beam epitaxy (MBE). In 
original SOS the system is represented as a simple 
cubic crystal. No vacancies and overhang are allowed. In 
this case, surface may be described as two-dimensional 
matrix. Indexes of the matrix serve as spatial coordi- 
nates of the substrate. Value of the matrix element is 
the height - number of atoms above the substrate in col- 
umn at this particular coordinate. Usually, three kind of 
processes are considered during the growth: deposition, 
surface diffusion and evaporation. We do not allow atoms 
to evaporate in our model. 

Deposition of atoms occurs at a random location on 
the substrate with the rate 
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(1) 



where F is the incoming beam fiux of deposited material 
and a is the distance between two nearest neighbors on 
the lattice. 

Surface diffusion is described by the set of hopping 
rates given by the Arrhenius distribution: 

-Ei/ksT 



(2) 
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where is a hopping rate of adatom in configuration z, 
vq - attempt frequency, fc^ is the Boltzmann's constant, 
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T - temperature of the substrate. Energy barrier Ei is 
calculated by the simple bond-counting scheme with the 
generic form: 

J 

E, = Es+Y,4^Ej, (3) 
i=i 

where Eg is the interaction energy of the free adatom 
with the substrate material. J is the total number of 
neighbors considered, n^*'' G (0,1), j = I,..., J is the 
number of bonds to the considered neighbor (which may 
be 0, when there is no neighbor, and 1 otherwise) in a 
configuration i. Ej are corresponding binding energies. 

In our model, we consider interactions between the 
adatom and its lateral first, second and third nearest 
neighbors. In the case of bcc(llO) surface, there are 4 
first nearest, 2 second nearest, and 2 third nearest neigh- 
bors (see Figure [1]). Hence, J — 8. Interaction energies 
Ej are material-related parameters of the model. Fur- 
ther in the text we will use notations £'2^, E^n for 
interaction energies with first, second and third nearest 
neighbors correspondingly. Energies for the Fe on Mo 
system were calculated by Yang and Asta [<!] employ- 
ing spin-polarized electronic density functional theory, 
generalized-gradient approximation 7\ calculations. De- 
tails of these calculations will be discussed in a future 
publication. Other model parameters do not depend on 
the material: substrate temperature, incoming flux of the 
deposited material and the surface coverage. 

In more detailed models, such effects as possibility of 
edge diffusion of an adatom, influence of a step-edge on 
interlayer transport, or presence of the strain (due to 
lattice mismatch between substrate and deposited ma- 
terial), should be taken into account. However, these 
effects are not consider at this stage of our work. 

B. Algorithm 

The most common way to simulate MBE growth is to 
use Kinetic Monte Carlo algorithm, also referenced to as 
BKL. It was proposed by Bortz et al. in 1975 [S]. In the 
core of KMC algorithms lie two key mechanisms. The 
first is the keeping of the list of all possible events in the 
system. The second one is the selecting particular event 
by making a search through the list of events. Since 
the appearance of BKL paper several different versions 
of KMC algorithm have been developed, for review see 
e.g. 9]. Different versions use different approach to or- 
ganizing event lists and different search algorithms. The 
most used KMC algorithms are the algorithm with linear 
search (LS), binary search (BS) (which is useful in case 
of large number of events and differs from LS algorithm 
by the method of bookkeeping: the list of events is orga- 
nized in form of the binary tree), K- level search [l^] and, 
finally, binning method by Maksym [lll |. 

For our simulations we adopt Maksym's method mod- 
ified by Haider et al. which allows to prevent events 
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FIG. 1: Fragment of the bcc(llO) surface with allowed direc- 
tions of jumps of the adatom (dark circle) . Lateral neighbors 
that interact with the adatom (nearest, second nearest and 
third nearest neighbors) are shown as light-gray circles and 
are marked by corresponding roman numerals. Arrows show 
allowed directions of the jump. 

with high hopping rates from dominating the simulation. 
In the following, we give short description of the algo- 
rithm with implemented specifications to suit our model. 

Fragment of bcc(llO) surface is shown on Figurc[T] As 
was mentioned before, the surface is represented as two- 
dimensional array with indices being the coordinates of 
the substrate site and the value of the array element is 
the number of deposited adatoms at that location. To 
eliminate any preferences in selecting particular sites we 
create a structure, that contains pointers to the random 
lattice sites. In order to speed up the search procedure, 
lattice is divided into a number of smaller groups. Each 
group contains pointer to the first element of the site list 
in this group, and number of events of a certain type. 
And again, to avoid any preferences in selecting certain 
group, we create a structure that contains pointers to 
random groups. 

In order to select particular event and location of this 
event, first we need to calculate total rate of all events in 
the system. We calculate this rate as a sum of all possible 
event rates weighted with the number of atoms in each 
local configuration type: 

M 

i?(A^i,...,7VM) = A^r-F^A^„r„, (4) 

a=l 

where N is the total number of surface atoms, r is the 
deposition rate ([T]), M is the total number of configura- 
tion types, Na - number of atoms with configuration type 
a, - hopping rate of the configuration a In our 
algorithm, individual hopping rates are calculated using 
the following diffusion barrier of the adatom: 

E = Es+ niEn + ^2,1 + nsEsn, (5) 

where ni, 712, are the counts of number of occupied 
first, second and third nearest neighbors correspondingly. 
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Taking into account geometry of the lattice (Figure [T]), 
there are 5 combinations of the nearest neighbors con- 
figurations, in which number of neighbors can be from 
to 4, and 3 combinations of second and third nearest 
neighbors, where number of neighbors can be from to 
2. Hence, total number of different configuration types 
M = 45. 

Individual probabilities of events may be calculated as 
follows: 



P(r) = Nt/R, 



(6) 



where P{t) is the probability of the deposition event and 
P{Ta) is the probability of the hopping event with config- 
uration a. Then a cumulative probability table is defined: 



C„ = ^P(r,), a = l,...,Af 
Cm+i^Cm+P{t). 



(7) 



In order to select an event a random number 5Ri in the 
region [0, R) is generated. Event a is then selected for 
which 



Cq-1 < 5Ri ^ Cq, for hopping event. 



Cm < 5Ri s$ C. 



M+l- 



for deposition event. 



(8) 



From this equation we can see, that if Sfti > Cm then 
deposition event will be realized. In other case one of the 
hopping events will be selected. 

The last step is to select the lattice site at which cho- 
sen event occurs. In case of the hopping event it is done 
in two substeps. First - search through the list of groups 
to find a group, that contains non-zero amount of events 
of a given type. To suppress any preferences in select- 
ing groups, every time the search starts from the random 
group. After the group is found, the search continues in- 
side the group (also, every time it starts from the random 
element of the group). The destination site for the jump 
is chosen randomly out of 6 adjacent neighbors. In our 
algorithm, we allow adatom to jump to its lateral first 
and second nearest neighbors (see Figure [T]). This choice 
was made mainly because of the geometrical considera- 
tions: first and second nearest neighbors are almost at 
the same distance from the jumping adatom. Besides, 
allowing atoms to jump to the third nearest neighbors 
sites would considerably increase the complexity of the 
algorithm realization. In case of the deposition event, 
site is chosen randomly. 

After the performing an event, it is needed to update 
the simulation timer. Since events are determined by a 
uniformly distributed random numbers, time interval be- 
tween two consecutive events t is derived from the Pois- 
son distribution of time intervals: 




FIG. 2: Change of island shape with temperature growth: a) 
T = 325iC, b) T = 3507^, c) T = 4007^, d) T = 425is:, e) 
T = 450if . Conditions of simulations are given in the text. 



where R is the total rate of all events. From this dis- 
tribution we can derive the time interval between events 



t 



log^ 
R 



(10) 



P{t)dt = Re^'^^dt 



(9) 



where 3^2 is a random number between and 1, i? - total 
rate of all possible events. 

Simulation continues until the terminal value of sur- 
face coverage is reached. Particular conditions of each 
simulation will be given during their discussion. 



III. RESULTS 

A. Dependence of island morphology on the 
substrate temperature 

In order to evaluate change of island shape, we need to 
define some quantitative parameter, therefore, we intro- 
duce average aspect ratio (a.a.r.). It is calculated as an 
average of ratio of transversal and longitudinal dimen- 
sions for each island averaged over all atomic islands in 
the system. Another quantity, that we will use in our 
discussion is the island compactness. It will be deter- 
mined visually - if island has a shape of a fractal then 
it is not compact. A fractal is a rough or fragmented 
geometric shape that can be subdivided in parts, each of 
which is, in a statistical sense, a reduced-size copy of the 
whole. Fractals are generally self-similar and indepen- 
dent of scale. We will refer to the fractals that are less 
fragmented as to fat fractals. We will call island compact 
if it is not a fractal. 
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e=0.01 ML 




Average Aspect Ratio vs. Temperature 



6=0.05 ML 



FIG. 3: Color online. Island shapes and evolution with cov- 
erage for two temperature values, a.) T — 500K, calculated 
energies; b) T = 3507^, scaled energies (see explanation in 
the text). Fragment of the lattice was scaled down. 



Figure [5] shows examples of the islands that were 
grown for different temperature of the substrate. Ex- 
cept for temperature, all other simulation parameters 
were the same: lattice size 300 x 300, incoming beam 
flux F — 0.01 ML/s (monolayers per second; this value 
is constant for all simulations), final surface coverage 
9 = 0.05 ML, atomic vibrational frequency Fmi, = 4 x 
10^^ Hz, and interaction energies are Eg = 0.4 eV, En — 
0.329 eV, E2n = 0.072 eV, Sgn = 0.079 eV. 

Hopping rates of the atoms increase with the temper- 
ature growth. Consequently, overall time spent for the 
surface diffusion also is growing with the temperature. 
That means that during the deposition on the substrate 
with the higher temperature atoms will have more time 
to diffuse to energetically preferable locations. Those lo- 
cations are the ones with the maximum number of lateral 
neighbors. Hence, for higher temperatures islands are be- 
coming more compact. 

Continuation of the shape evolution is shown on the 
Figure [31 First island snapshot is the island at T = 50QK 
and the same conditions described above. To eliminate 
finite size effect, lattice size was increased to 1000 x 1000. 
We can see that an island is clearly bordered by two kinds 
of facets. Figure shows the island obtained under an- 
other conditions. For this simulation all energies and the 
substrate temperature were scaled by a factor a = 0.5. 
Initially, scaling was done in order to reduce computa- 
tion time. Multiplying both energies and temperature 
by the same factor does not change the hopping rates of 
the atoms Hence, we may say that the same result 
would be obtained for initial set of energies (Figures[2]and 
[3^) and temperature T = TOOK. Figure [3] was scaled to 
fit the picture. Originally, the size of this rounded island 
is much larger that the size of rhombic-like island on the 
Figure [3^. 

One can observe three stages of the shape change on 
the FiguresmandU from T = 325K to T = AOOK islands 
are changing from being fractals to fat fractals. After 
that the islands shape is becoming compact quadrangle 
at T = 500-ftr. As temperature growth further, islands 
start to loose their distinc facets and obtain rounded 




Insets: examples of Island shape 
at the particular temperatures 
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FIG. 4: Dependence of the average aspect ratio on the sub- 
strate temperature. Insets show examples of island shapes 
at different temperatures. Conditions of the simulations are 
given in the text. 



shape as shown on the Figure ^jp. Similar islands with 
rounded shapes were observed in the experimental work 

MM- 

Another type of shape evolution of the islands may be 
observed on the Figure [31 evolution of the island shape 
with the coverage. We can observe, that islands obtain 
their shape at the very early stages of the growth. As 
the deposition proceeds, islands only grow in size and do 
not change the shape. Of course, there is a possibility for 
two or more islands to merge. In this case (shortly after 
the merge) the resulting island will not have the shape of 
other individual islands. But given enough time for the 
atoms to diffuse, merged island shape will evolve into 
rhombic-like or rounded island. This is not the case for 
the fractal and fat fractal islands. Since atoms are less 
mobile on the substrate with low temperatures, fractal 
islands will remain fractals after the merge. 

Figure [4] shows dependence of the a.a.r. on the sub- 
strate temperature. Data for this graph was obtained 
under the following simulation conditions: for temper- 
atures T 180..200ii: lattice size was 300 x 300, for 
T = 210..230ii: - 600 x 600 and for T = 2A0K size was 
1000 X 1000, final surface coverage 6 = 0.1 AIL; atomic 
vibrational frequency Fyn, = 2 x 10^^ Hz; interaction 
energies Es = 0.2 eV, E^ = 0.14 eV, E2n = 0, E^^ = 
0.06 eV. For 300 x 300 lattices data was averaged over 
10 independent simulations and for larger lattices - over 
5 independent runs. 

Dependence of average aspect ration on the tempera- 
ture that is shown on the Figure [4] supports our previ- 
ous observations on the island shape. As temperature 
growth, islands' average aspect ratio is growing, indicat- 
ing that islands becoming more compact. This statement 
is also illustrated by the insets on the Figure [H These 
insets show examples of the islands at different tempera- 
tures and we can see the transition from the islands with- 
out sharp facets to the compact rhombic-like islands. 



5 



B. Dependence of island morphology on the 
interaction energies 

Temperature is not the only parameter, that has an 
influence on the island shape. In this section we present 
the results on how the interaction energies of the adatom 
with its lateral neighbors change the islands shape. Since 
in the following simulations we are changing material- 
dependent parameters (interaction energies), this section 
is not strictly related to the Fe/Mo(110) growth process. 
Nevertheless, these results may give usefuU insight to the 
island morphology for the systems with other sets of in- 
teraction energies. 

Average Aspect Ratio vs.Ejn 



Average Aspect Ratio vs. Ejn 



Insets: examples of the island shape 
at the particular values of E2n 





(E2n^0.14,eV) 



0.04 
E2„ [=V] 



FIG. 5: Dependence of average aspect ratio on interaction 
energy with second nearest neighbor. Insets show islands at 
different values of E2n 



1 

0.9 

D, 



0, 
0, 
n 

D, 
0, 
0, 
0, 



Insets: examples of the Island shape at the 
particular values of Ean 
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FIG. 6: Average aspect ratio as a function of interaction en- 
ergy with third nearest neighbor. Insets show islands at dif- 
ferent values of Ezn 



stronger than that of E2n ■ Although one can observe the 
decrease of a.a.r. with growth of E^n (as in the case of 
changing i?2n), the change of island shape differs from the 
case above. Here with the change of E^n islands become 
more prolongated which leads to the decrease of average 
aspect ratio. But at the same time the compactness of 
the islands is preserved - they are not becoming fractals 
or even fat fractals in the range of simulations. Change 
of the island shape is also illustrated by the insets on the 
Figure O where one can see islands at different values of 



Figureinishows the dependence of average aspect ration 
on the interaction energy of the adatom with the second 
neareast neighbors. To obtain this result, the following 
simulation conditions were used: 300 x 300 lattice, 9 = 
0.1 ML, T = 230X, F„6 = 2 x lO^^ jj^^ interaction 
energies: Es = 0.2 eV, E^ = 0.14 eV, Es^ = 0.03 eV. 
The data was averaged over 10 independent runs. 

As one can see, the main influence of i?2n is not the 
change of a.a.r., but the change of shape of the islands. 
The shape is changing from compact at low values of i?2n 
to the shape without sharp facets and, eventually, to the 
fractals. This evolution is illustrated by the insets on 
the Figure [5] which are showing snapshots of the islands 
at different values of i?2n- Simulations were terminated 
at E2n — 0.08 eV since around that value island shape 
started to become fat fractal and measuring a.a.r. is 
usefuU only for compact islands. 

Dependency of the island morphology on E^n is shown 
on the Figure [HI Here we used 300 x 300 lattice, terminal 
coverage 9 — 0.1 ML, substrate temperature T — 250K, 
atomic vibrational frequency Fyn, = 2 x 10^^ Hz and in- 
teraction energies Es = 0.2 eV, En = 0.14 eV, i?2n = 
0.04 eV. Each data point was averaged over 10 indepen- 
dent runs. 

Influence of E^n on the island shape appears to be 



IV. CONCLUSIONS 

In this paper, we have presented results of KMC sim- 
ulations of sub-monolayer island growth at bcc(llO) sur- 
face. We have utilized extended bond counting model 
with interactions up to third nearest neighbors. Interac- 
tions were based on ab initio calculations. However, for 
this parameter set we were able to perform simulations 
only for temperatures up to bOOK. In order to obtain re- 
sults for higher temperatures we used scaled interaction 
energies. 

We observed usual transition of island shape from frac- 
tal to compact. Compact islands at T = 5QQK have 
rhombic-like shape. For higher temperatures they obtain 
rounded shape in agreement with the experiment p^lT5|. 
In order to quantify morphological change we have eval- 
uated the dependence of average aspect ratio of island 
on the temperature and values of adatom interactions. 
We observed that with increasing temperature aspect ra- 
tio is increasing, and eventually saturates. In the case 
of dependence on the strength of interaction, we found 
that the aspect ratio decreases slowly with the increase 
of interaction energy with the second neighbor. But the 
island shape is quickly transforming from compact to fat 
fractals. With the increase of interaction energy with the 
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third neighbor, the aspect ratio is decreasing much faster 
than in the previous case. At the same time, the island 
shape remains compact in the range of simulation. 

In addition, we observed that the system may have a 
state where the islands develop an extra facet and re- 
semble hexagon-like shapes. This our result is similar to 
the experimental works where hexagon islands were ob- 
served [l6i] . However, the conditions (energies, temper- 
atures and coverage) under which islands with hexagon- 
like shapes are grown are still needed to be clarified. In 
the case of heteroepitaxy, it is also important to con- 
sider effect of strain which is currently not included in 
our model. In order to obtain more realistic results we 



are working on the new algorithm which will include edge 
diffusion and strain. 

We have realized that our calculations for high tem- 
peratures (T > 500K) and interaction energies obtained 
by ab initio calculation are very demanding in terms of 
computer time. Nevertheless, our results are in qualita- 
tive agreement with experimental findings. 
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